/* f_T8_row_joint_pvals.do *****************************************************

For Table 8 (Subgroup analysis), find the p-values for a test
where the null hypothesis is that there are no differences in estimates across
all subgroups for a given outcome

*******************************************************************************/

* Note: Since ivreghdfe is incompatible with suest, we estimated SUR by hand,
* by stacking the dataset

estimates clear

* Switches
local non_mcas 		1
local mcas			  1

include "${code}/ado/personal/group_p_score.ado"

* Section I: Do all non-MCAS outcomes
if `non_mcas' == 1 {
	use "$stata_data_analysis/final_sample.dta", clear
	drop grade

	gen male = female == 0

	local outcome_list any_6m_asgn any4yr_nsc_6m_asgn  ///
					   everany_nsc everany4yr_nsc  has_sat_took_10 ///
					   grad_hs behav_index_tot ///


	local group_vars "frpl not_frpl male female black white hispanic"


	matrix pval = J(60,6,.)
	local i = 1

	gen id = _n

	* Loop over subgroups
	local subgroup_count = 1
	expand 7
	gen y = .
	bysort id: gen model = _n

	* Append MCAS stacked dataset:

	forvalues y = 1/7{
		gen model`y' = 1 if model == `y'
		display in red "subgroup: `y'"


			display in red "group var: `group'"
			foreach v of varlist att soffany pformany {
				gen `v'_model`y' = `v' * (model == `y')
			}

	}

	drop if model==1 & frpl!=1 // frpl
	drop if model==2 & not_frpl!=1 // not frpl
	drop if model==3 & male!=1 // male
	drop if model==4 & female!=1 // female
	drop if model==5 & black!= 1 // black
	drop if model==6 & hispanic!=1 // hispanic
	drop if model==7 & white!=1	// white

	foreach y of local outcome_list {
		* Joint regression
		ivreghdfe `y' (att_model* = soffany_model*) if in_K1 & !inlist(pformany, 0, 1), robust absorb(pformany_model*) // dummy version
		estimates store joint`y'

		qui: cap estadd scalar N = r(N)

		test (_b[att_model1] = _b[ att_model2])	///
			 (_b[att_model3] = _b[ att_model4]) ///
			 (_b[att_model5] = _b[ att_model6] = _b[att_model7])

		qui: cap estadd scalar f = r(F)
		qui: cap estadd scalar p = r(p)

		* Export results
		cd "$results/tables/joint_f_outcome/"
		print_results  att_model* , model_names( joint`y' ) hidevars_print( 0 ) outfile("joint`y'.csv") ///
			stat_names("N" "f" "p") stat_labs("Sample Size" "F test" "p-value" )
	}
}

* Section II: Do all non-MCAS outcomes
if `mcas' == 1 {
	use "$stata_data_analysis/full_samp.dta", clear

	keep studentno att in_sims in_K1 soffany yr_asgn hispanic asian female black white ///
	age_PreK_enr num_ranked first_walk pop med_inc_fam pov share_black share_white share_hispanic ///
	std_scaled_e3-std_scaled_e8 ///
	std_scaled_e10 std_scaled_m3-std_scaled_m8 std_scaled_m10 ELC soff* pform* ///
	frpl not_frpl nLC sp_biling inc1 inc2 inc3 inc4  ///
	std_e3-std_e8 std_e10 std_m3-std_m8 std_m10

	gen male = female == 0

	reshape long std_scaled_e std_scaled_m behav_index, i(studentno) j(grade)

		local i = 1

	gen id = _n

	* Loop over outcomes
	local subgroup_count = 1
	expand 7
	gen y = .
	bysort id: gen model = _n

	forvalues y = 1/7{
		gen model`y' = 1 if model == `y'
		display in red "subgroup: `y'"


			display in red "group var: `group'"
			foreach v of varlist att soffany pformany {
				gen `v'_model`y' = `v' * (model == `y')
			}

	}

	drop if model==1 & frpl!=1 // frpl
	drop if model==2 & not_frpl!=1 // not frpl
	drop if model==3 & male!=1 // male
	drop if model==4 & female!=1 // female
	drop if model==5 & black!= 1 // black
	drop if model==6 & hispanic!=1 // hispanic
	drop if model==7 & white!=1	// white

	local outcome_list std_scaled_m

	foreach y of local outcome_list {
		* Joint regression
		ivreghdfe `y' (att_model* = soffany_model*) if in_K1 & !inlist(pformany, 0, 1), robust absorb(pformany_model*) cluster(studentno)  // dummy p-score

		estimates store joint`y'

		qui: cap estadd scalar N = r(N)

		test (_b[att_model1] = _b[ att_model2])	///
			 (_b[att_model3] = _b[ att_model4]) ///
			 (_b[att_model5] = _b[ att_model6] = _b[att_model7])

		qui: cap estadd scalar f = r(F)
		qui: cap estadd scalar p = r(p)

		* Export results
		cd "$results/tables/joint_f_outcome/"
		print_results  att_model* , model_names( joint`y' ) hidevars_print( 0 ) outfile("joint`y'.csv") ///
			stat_names("N" "f" "p") stat_labs( "F test" "p-value" )
	}
}
